
clc 
clear
fsize = 14;
format short g                                      % display format for all numbers in consolve
set(0,'DefaultTextInterpreter','latex');            % set text interpreter for figures
set(0,'DefaultLegendInterpreter','latex');          % set text interpreter for all legends in figures
set(0,'DefaultAxesTickLabelInterpreter','latex');   % set text interpreter for axes in figures
set(0,'DefaultAxesFontSize',fsize);                 % set font size for axes in figures
set(0,'DefaultTextFontsize',fsize);                 % set font size for other text
set(0,'DefaultLegendFontsizeMode','manual');        % set font size for legend
set(0,'DefaultLegendFontsize',fsize);               % set font size for other text
set(0,'DefaultLineLineWidth',2)                     % set line thickness

green = [0,0.5,0];
white = [1,1,1];
grey = 0.6 * white;


%% file paths and names for read in and saving
RootPath = 'C:\Users\dempsey.164\Dropbox\Lending_standards_and_credit\Final_Replication\';
DataPath = strcat(RootPath, 'Data\');
FigPath = strcat(RootPath, 'Figures\');
ModelPath = strcat(RootPath, 'FortranOutput\');

BasePath = strcat(ModelPath, 'Baseline\');
ShortTermPath = strcat(ModelPath, 'ShortTerm\');
NoIntCostPath = strcat(ModelPath, 'LowIntCost\');
FixedRecoveryPath = strcat(ModelPath, 'BKOnly\');

Paths = {BasePath, ShortTermPath, NoIntCostPath, FixedRecoveryPath};
n = length(Paths);

PanelPath = strcat(BasePath,'Panel\');

%% load in data
% balance metrics for Figures 1 and 3
XDataFile = strcat(DataPath,'dataforfigures_Sept2024.csv');
X = readtable(XDataFile);
X = X(1:end-1,:);           % eliminate the last bucket which is a catch-all outlier
Bins = X.even_fine;
ActualDefaultProb = 100 * X.mean_default_actual_2;
PredictedDefaultProb = 100 * X.mean_default_predict_2;
Count = X.count;
ShareBalances_Revolved = X.share_revolved_balance;
ShareBalances = X.share_balance;
AvgBalance = X.avg_revolved_balance;
ShareRevolvers = 100 * X.revolver_share;
ShareRevolversQ1 = 100 * X.q1_revolver;
ShareRevolversQ2 = 100 * X.q2_revolver;
ShareRevolversQ3 = 100 * X.q3_revolver;
ShareRevolversQ4 = 100 * X.q4_revolver;

% rescale to the right level of default
MaxDefault = 20;
i = gridlookup(length(Bins), Bins, MaxDefault);
ActualDefaultProb = ActualDefaultProb(1:i);
PredictedDefaultProb = PredictedDefaultProb(1:i);
Count = Count(1:i);
ShareAccounts = 100 * Count / sum(Count);
ShareBalances = 100 * ShareBalances(1:i) / sum(ShareBalances(1:i));
ShareBalances_Revolved = 100 * ShareBalances_Revolved(1:i) / sum(ShareBalances_Revolved(1:i));
AvgBalance = AvgBalance(1:i);
ShareRevolvers = ShareRevolvers(1:i);
ShareRevolversQ1 = ShareRevolversQ1(1:i);
ShareRevolversQ2 = ShareRevolversQ2(1:i);
ShareRevolversQ3 = ShareRevolversQ3(1:i);
ShareRevolversQ4 = ShareRevolversQ4(1:i);




% spread metrics for Figure 2 
XDataFile = strcat(DataPath,'stats_even_groups_full_fine_revolved.csv');

X = readtable(XDataFile);
X = X(1:end-1,:);           % eliminate the last bucket which is a catch-all outlier

MeanAPR = X.mean_apr;
SDAPR = X.sd_apr;
MeanSpread = 100 *  X.mean_spread_a;    % multiply by 100 to convert to percentage points
SDSpread = 100 * X.sd_spread;

% rescale to the right level of default
MeanAPR = MeanAPR(1:i);
SDAPR = SDAPR(1:i);
MeanSpread = MeanSpread(1:i);
SDSpread = SDSpread(1:i);




%% load in common model features 
kappa = f2m(strcat(BasePath,'kappa'),1);
phi = f2m(strcat(BasePath,'phi'),1);

% grid sizes
na = f2m(strcat(BasePath,'na'),1);
naNeg = f2m(strcat(BasePath,'naNeg'),1);
aZeroInd = naNeg + 1;
ni = f2m(strcat(BasePath,'ni'),1);
ne1 = f2m(strcat(BasePath,'ne1'),1);
ne2 = f2m(strcat(BasePath,'ne2'),1);
ne3 = f2m(strcat(BasePath,'ne3'),1);

wage = f2m(strcat(BasePath,'wage'),1);

% grids
aGrid = f2m(strcat(BasePath,'aGrid'),na);
e1Grid = f2m(strcat(BasePath,'e1Grid'),ne1);
e2Grid = f2m(strcat(BasePath,'e2Grid'),ni);
e3Grid = f2m(strcat(BasePath,'e3Grid'),ni);
yMat = zeros(ni, ne1);
for e1i = 1:ne1
for i = 1:ni
   yMat(i,e1i) = wage * e1Grid(e1i) * e2Grid(i) * e3Grid(i);
end
end






%% Figure 1: Default and balances 
FigName = strcat(FigPath,'Fig_1_DefaultDistribution');

nr = 1;
nc = 3;
f = figure('Position',[0,0,nc * 450,nr * 300]);
xl = [0,MaxDefault];


subplot(nr,nc,1)
plot(PredictedDefaultProb, ActualDefaultProb,'k-o');
hold on 
plot(PredictedDefaultProb, PredictedDefaultProb,':','Color',grey);
xlim(xl)
title('(a) actual vs predicted default')
ylabel('realized default rate (\%)')
xlabel('predicted default probability (\%)')
hleg = legend('realized default rate',...
    '$45^\circ$ (pred. = realized)');
set(hleg, 'location','best','box','off')

% panel (b): share of accounts
subplot(nr,nc,2)
plot(PredictedDefaultProb, cumsum(ShareAccounts),'k-o')
hold on
plot(PredictedDefaultProb, cumsum(ShareBalances),'b-s')
plot(PredictedDefaultProb, cumsum(ShareBalances_Revolved),'r-.d')
ylabel('cumulative share (\%)')
xlabel('predicted default probability (\%)')
title('(b) distribution of accounts and balances')
hleg = legend('accounts','cycle-end bals.','revolved bals.');
set(hleg, 'location','southeast','box','off')

% panel (c): share revolving
subplot(nr,nc,3)
plot(PredictedDefaultProb, ShareRevolversQ1,'k-s')
hold on
plot(PredictedDefaultProb, ShareRevolversQ2,'b-.d')
plot(PredictedDefaultProb, ShareRevolversQ3,'r:o')
plot(PredictedDefaultProb, ShareRevolversQ4,'Color',green,'Linestyle','--','Marker','x');
xlim(xl)
title('(c) share of accounts revolving')
xlabel('predicted default probability (\%)')
ylabel('share (\%)')
hleg = legend('at least 1 qtr', 'at least 2 qtrs','at least 3 qtrs','at least 4 qtrs');
set(hleg, 'location','best','box','off')

saveas(f,strcat(FigName,'.eps'),'epsc')


%% Figure 2: Spreads
FigName = strcat(FigPath,'Fig_2_Spreads');

% linear fit to the data
one = ones(length(PredictedDefaultProb),1);
mdl = fitlm(PredictedDefaultProb,MeanSpread);
intercept = mdl.Coefficients.Estimate(1);
intercept_SE = mdl.Coefficients.SE(1);
slope = mdl.Coefficients.Estimate(2);
slope_SE = mdl.Coefficients.SE(2);
AdjustedR2 = mdl.Rsquared.Adjusted;

% theoretical prediction
rf = 0.055;     % benchmark / risk-free rate 
iota = 0.05;    % intermediation cost
Theoretical_ST_NoRecovery = 100 * ((1.0+rf+iota) * one ./ (one - PredictedDefaultProb / 100) - (1.0+rf) );
Theoretical_ST_FullRecovery = 100 * iota * one;

% FIGURE
f = figure('Position',[0,0,1200,400]);
nr = 1;
nc = 2;
xl = [0,MaxDefault];


subplot(nr,nc,1)
a = plot(PredictedDefaultProb, MeanSpread,'k-o');%,'Color',col1);
hold on
b = plot(ActualDefaultProb, MeanSpread,'b:');%,'Color',col2);
c = plot(PredictedDefaultProb, intercept + slope * PredictedDefaultProb,'r-.');%,'Color',col3);
xlim(xl)
text(13.5,17.2,{...
    ['intercept ',num2str(round(intercept,2)),' (',num2str(round(intercept_SE,2)),')'],...
    ['slope ',num2str(round(slope,2)),' (',num2str(round(slope_SE,2)),')']},...
    'Color','red','FontSize',10)
title('(a) spreads vs predicted default')
ylabel('average spread (annualized pp)')
xlabel('default probability (\%)')
hleg = legend([a,b,c],'average spread vs. predicted',...
    'average spread vs. realized',...
    ['linear fit (adj. $R^2$ = ',num2str(round(AdjustedR2,3)),')']);
set(hleg, 'location','southeast','box','off')


subplot(nr,nc,2)
basevalue = min(Theoretical_ST_FullRecovery);
d = area(PredictedDefaultProb, Theoretical_ST_NoRecovery,'FaceColor','b','EdgeColor','b','FaceAlpha',.1,'Linewidth',.01,'basevalue',basevalue);
hold on
area(PredictedDefaultProb, Theoretical_ST_FullRecovery,'FaceColor',white,'EdgeColor','b','FaceAlpha',1,'Linewidth',.01,'basevalue',basevalue);
b = plot(PredictedDefaultProb, Theoretical_ST_NoRecovery,'r-');
c = plot(PredictedDefaultProb, Theoretical_ST_FullRecovery,'-','Color',green);
a = plot(PredictedDefaultProb, MeanSpread,'k-o');%,'Color',col1);d = plot(PredictedDefaultProb, Theoretical_ST_NoRecovery,'--','Color',col3);%,'Color',col3);
xlim(xl)
yl = ylim;
ylim([0,yl(2)])
text(.5,100*iota*0.7,'$\nwarrow$ intermediation cost')
title('(b) comparison with short term debt model')
ylabel('average spread (annualized pp)')
xlabel('default probability (\%)')
hleg = legend([a,b,c,d],'data',...
    'theory (short term, no recovery)',...
    'theory (short term, full recovery)',...
    'range spanned by recovery');
set(hleg, 'location','northwest','box','off')

saveas(f,strcat(FigName,'.eps'),'epsc')


%% Figure 3: Revolving behavior, model v data 
FigName = strcat(FigPath,'Fig_3_Revolvers_ModelVData');

% data profile
CondShare_Q2 = ShareRevolversQ2;
CondShare_Q3 = ShareRevolversQ3;
CondShare_Q4 = ShareRevolversQ4;

DataReg = fitlm(PredictedDefaultProb, CondShare_Q2);
DataQ2Coefs = DataReg.Coefficients.Estimate;
DataQ2Fit = DataQ2Coefs(1) + DataQ2Coefs(2) * PredictedDefaultProb;

DataReg = fitlm(PredictedDefaultProb, CondShare_Q3);
DataQ3Coefs = DataReg.Coefficients.Estimate;
DataQ3Fit = DataQ3Coefs(1) + DataQ3Coefs(2) * PredictedDefaultProb;

DataReg = fitlm(PredictedDefaultProb, CondShare_Q4);
DataQ4Coefs = DataReg.Coefficients.Estimate;
DataQ4Fit = DataQ4Coefs(1) + DataQ4Coefs(2) * PredictedDefaultProb;


% model profile 
nSim = f2m(strcat(PanelPath,'nSim'), 1);
tSim = f2m(strcat(PanelPath,'tSim'), 1);
nObs = nSim * tSim;

DefProb = f2m(strcat(PanelPath,'DefProbSim'), nSim);

x = [nSim,tSim];
Assets = f2m(strcat(PanelPath,'aSim'), x);

nPartition = 40;
pDPartition = linspace(0.0,0.2, nPartition+1)';

% compute revolver status in the model
nq = 4;
RevolverStatus = zeros(nSim, tSim-1, nq, nPartition);
RevolverShare = zeros(nq,nPartition);
for j = 1:nPartition
    BinFlag =  (DefProb >= pDPartition(j)) .* (DefProb < pDPartition(j+1));
    nCount = sum(BinFlag);
    for i = 1:nq
        for t = 2:tSim
            if (t <= tSim + 1 - i)
                RevolverStatus(:,t,i,j) = (sum(Assets(:,t:t+i-1) <= 0,2) >= i).* BinFlag;
            end
        end
        RevolverShare(i, j) = 100 * sum(RevolverStatus(:,:,i,j),[1,2]) / (nCount * (tSim - i));
    end
end

Model_CondShare_Q2 = 100 * RevolverShare(2,:) ./ RevolverShare(1,:);
Model_CondShare_Q3 = 100 * RevolverShare(3,:) ./ RevolverShare(1,:);
Model_CondShare_Q4 = 100 * RevolverShare(4,:) ./ RevolverShare(1,:);

pDPlot = 100.0 * pDPartition(2:end);

ModelReg = fitlm(pDPlot, Model_CondShare_Q2);
ModelQ2Coefs = ModelReg.Coefficients.Estimate;
ModelQ2Fit = ModelQ2Coefs(1) + ModelQ2Coefs(2) * pDPlot;

ModelReg = fitlm(pDPlot, Model_CondShare_Q3);
ModelQ3Coefs = ModelReg.Coefficients.Estimate;
ModelQ3Fit = ModelQ3Coefs(1) + ModelQ3Coefs(2) * pDPlot;

ModelReg = fitlm(pDPlot, Model_CondShare_Q4);
ModelQ4Coefs = ModelReg.Coefficients.Estimate;
ModelQ4Fit = ModelQ4Coefs(1) + ModelQ4Coefs(2) * pDPlot;

% figure
yl = [20,100];
nr = 1;
nc = 3;
sz = 25;
f = figure('Position',[0, 0, nc * 450, nr * 300]);

subplot(nr,nc,1)
scatter(PredictedDefaultProb, CondShare_Q2,sz,'o','MarkerEdgeColor','k')
hold on
scatter(pDPlot, Model_CondShare_Q2,sz,'d','MarkerEdgeColor','b')
plot(PredictedDefaultProb, DataQ2Fit, 'k--') 
plot(pDPlot, ModelQ2Fit,'b--')
xlim(xl)
ylim(yl)
title('(a) share $\ge$ 2Q given 1Q')
xlabel('default probability (\%)')
ylabel('\%')
hleg = legend(['data (int = ',num2str(round(DataQ2Coefs(1),1)),', slope = ',num2str(round(DataQ2Coefs(2),1)),')'],...
    ['model (int = ',num2str(round(ModelQ2Coefs(1),1)),', slope = ',num2str(round(ModelQ2Coefs(2),1)),')']);
set(hleg, 'location','best','box','off')


subplot(nr,nc,2)
scatter(PredictedDefaultProb, CondShare_Q3,sz,'o','MarkerEdgeColor','k')
hold on
scatter(pDPlot, Model_CondShare_Q3,sz,'d','MarkerEdgeColor','b')
plot(PredictedDefaultProb, DataQ3Fit, 'k--') 
plot(pDPlot, ModelQ3Fit,'b--')
xlim(xl)
ylim(yl)
title('(b) share $\ge$ 3Q given 1Q')
xlabel('default probability (\%)')
hleg = legend(['data (int = ',num2str(round(DataQ3Coefs(1),1)),', slope = ',num2str(round(DataQ3Coefs(2),1)),')'],...
    ['model (int = ',num2str(round(ModelQ3Coefs(1),1)),', slope = ',num2str(round(ModelQ3Coefs(2),1)),')']);
set(hleg, 'location','best','box','off')

subplot(nr,nc,3)
scatter(PredictedDefaultProb, CondShare_Q4,sz,'o','MarkerEdgeColor','k')
hold on
scatter(pDPlot, Model_CondShare_Q4,sz,'d','MarkerEdgeColor','b')
plot(PredictedDefaultProb, DataQ4Fit, 'k--') 
plot(pDPlot, ModelQ4Fit,'b--')
xlim(xl)
ylim(yl)
title('(c) share $\ge$ 4Q given 1Q')
xlabel('default probability (\%)')
hleg = legend(['data (int = ',num2str(round(DataQ4Coefs(1),1)),', slope = ',num2str(round(DataQ4Coefs(2),1)),')'],...
    ['model (int = ',num2str(round(ModelQ4Coefs(1),1)),', slope = ',num2str(round(ModelQ4Coefs(2),1)),')']);
set(hleg, 'location','best','box','off')

saveas(f,strcat(FigName,'.eps'),'epsc')




%% Figure 4: Spread function illustrations
FigName = strcat(FigPath,'Fig_4_SpreadIllustration');

% price schedules and distributions
spreadsAnnualized = zeros(naNeg, na, ni, ne1, n);
mu0 = zeros(na, ni, ne1, n);
for i = 1:n
    % load in price and distribution
    q = f2m(strcat(Paths{i},'q'), [naNeg, na, ni, ne1]);
    spreads = kappa * (q .^ (-1.0)) - phi;
    spreadsAnnualized(:,:,:,:,i) = 100 * ((1.0 + spreads) .^ 4.0 - 1.0);
    mu0(:,:,:,i) = f2m(strcat(Paths{i},'mu0'), [na, ni, ne1]);
end

% figure
nr = 1;
nc = 2;

ai = aZeroInd;

i = floor(ni/2);
i2 = i - 9;

e1i = 4;
e1i2 = e1i - 2;

f = figure('Position',[0,0,nc * 600, nr * 400]);
subplot(nr,nc,1)
plot(aGrid(1:naNeg), spreadsAnnualized(:,ai,i,e1i,1))
hold on
plot(aGrid(1:naNeg), spreadsAnnualized(:,ai,i2,e1i,1),':')
plot(aGrid(1:naNeg), spreadsAnnualized(:,ai,i,e1i2,1),'-.')
xlim([-0.5,0])
title('(a) varying income')
ylabel('spread (ann. pp)')
xlabel('debt choice $a^\prime$')
hleg = legend(['$\epsilon$ = (',num2str(round(e1Grid(e1i),1)),', ',num2str(round(e2Grid(i),1)),', ',num2str(round(e3Grid(i),1)),')'],...
    ['$\epsilon$ = (',num2str(round(e1Grid(e1i),1)),', ',num2str(round(e2Grid(i2),1)),', ',num2str(round(e3Grid(i2),1)),')'],...
    ['$\epsilon$ = (',num2str(round(e1Grid(e1i2),1)),', ',num2str(round(e2Grid(i),1)),', ',num2str(round(e3Grid(i),1)),')']);
set(hleg,'location','best','box','off')

subplot(nr, nc, 2)
plot(aGrid(1:naNeg), spreadsAnnualized(:,ai,i,e1i,1),'b-')
hold on
plot(aGrid(1:naNeg), spreadsAnnualized(:,ai,i,e1i,2),'Color',green,'Linestyle',':')
plot(aGrid(1:naNeg), spreadsAnnualized(:,ai,i,e1i,3),'r-.')
plot(aGrid(1:naNeg), spreadsAnnualized(:,ai,i,e1i,4),'Color',grey,'Linestyle','--')
xlim([-0.5,0])
title('(b) model variants')
ylabel('spread (ann. pp)')
xlabel('debt choice $a^\prime$')
hleg = legend('baseline','short term','low int. cost','bankruptcy only');
set(hleg,'location','south','box','on')

saveas(f,strcat(FigName,'.eps'),'epsc')






%% Table 7: Regression Validation

% data regressions come from Table 4 -- see that part of replication package
% load in model
nSim = f2m(strcat(PanelPath,'nSim'), 1);
tSim = f2m(strcat(PanelPath,'tSim'), 1);

DefProb = f2m(strcat(PanelPath,'DefProbSim'), nSim);

x = [nSim,tSim];
Assets = f2m(strcat(PanelPath,'aSim'), x);
Earnings = f2m(strcat(PanelPath,'ySim'), x);
Spreads = f2m(strcat(PanelPath,'SpreadSim'), x);
Spreads = 100 * ((1.0 + Spreads) .^ 4 - 1); % make annualized percentage terms 


% run regressions
tMax = 4;   % limit to 1 year sample
LowMidThresh = 0.005;
MidHighThresh = 0.1;
DefProbMat = repmat(DefProb, [1, tMax]);
IncomeAtOrig = repmat(Earnings(:,1), [1, tMax]);

% aggregate 
i = find((Assets(:,2:tMax+1) < 0) .* (Spreads(:,1:tMax) > 0));   % linear indices where borrowing happens;
DefProb_Reg = log(100*DefProbMat(i));
IncomeAtOrig_Reg = log(IncomeAtOrig(i));
Spreads_Reg = log(Spreads(i));

AggModel_1 = fitlm(DefProb_Reg, Spreads_Reg);
AggModel_2 = fitlm([DefProb_Reg, IncomeAtOrig_Reg], Spreads_Reg);

X = zeros(5,2);
X(1,1) = AggModel_1.Coefficients.Estimate(2);
X(2,1) = AggModel_1.Coefficients.SE(2);
X(5,1) = AggModel_1.Rsquared.Adjusted;

X(1,2) = AggModel_2.Coefficients.Estimate(2);
X(2,2) = AggModel_2.Coefficients.SE(2);
X(3,2) = AggModel_2.Coefficients.Estimate(3);
X(4,2) = AggModel_2.Coefficients.SE(3);
X(5,2) = AggModel_2.Rsquared.Adjusted;

writematrix(X,strcat(RootPath,'DI_JPEM_Model_Tables.xlsx'),'sheet','Table 7 Regression Validation','Range','C6:D10')







%% Figure 5: Spreads by default probability, main figure 
FigName = strcat(FigPath,'Fig_5_Spreads_ModelVData');
DataFile = strcat(DataPath,'stats_even_groups_full_fine_revolved.csv');

X = readtable(DataFile);
X = X(1:end-1,:);           % eliminate the last bucket which is a catch-all outlier

Bins = X.even_fine;
PredictedDefaultProb = 100 * X.mean_default_predict_2;
MeanSpread = 100 *  X.mean_spread_a;    % multiply by 100 to convert to percentage points

% rescale to the right level of default
MaxDefault = 20;
i = gridlookup(length(Bins), Bins, MaxDefault);
PredictedDefaultProb = PredictedDefaultProb(1:i);
MeanSpread = MeanSpread(1:i);

% linear fit
mdl = fitlm(PredictedDefaultProb,MeanSpread);
intercept = mdl.Coefficients.Estimate(1);
slope = mdl.Coefficients.Estimate(2);
DataLinearFit = intercept + slope * PredictedDefaultProb;

% load in model 
nPartition = 40;
ModelMeanSpread = zeros(nPartition, n);
ModelLinearFit = ModelMeanSpread;

ModelIntercept = zeros(n,1);
ModelIntercept_SE = ModelIntercept;
ModelSlope = ModelIntercept;
ModelSlope_SE = ModelIntercept;
ModelAdjustedR2 = ModelIntercept;
    
for i = 1:n
    % load in data 
    Spreads = f2m(strcat(Paths{i},'AvgBinSpread'), nPartition);
    if (i == 1)
        pDPartition = f2m(strcat(Paths{i},'pDPartition'), nPartition);
        pDPartition = pDPartition - 0.5;
    end
    Spreads(Spreads == 0) = NaN;
    ModelMeanSpread(:,i) = Spreads;
    
    % linear fit 
    mdl = fitlm(pDPartition,ModelMeanSpread(:,i));
    ModelIntercept(i) = mdl.Coefficients.Estimate(1);
    ModelIntercept_SE(i) = mdl.Coefficients.SE(1);
    ModelSlope(i) = mdl.Coefficients.Estimate(2);
    ModelSlope_SE(i) = mdl.Coefficients.SE(2);
    ModelAdjustedR2(i) = mdl.Rsquared.Adjusted;
    
    ModelLinearFit(:,i) = ModelIntercept(i) + ModelSlope(i) * pDPartition;
end



% FIGURE
xl = [0,MaxDefault];
nr = 1;
nc = 2;
figWidth = 600;
figHeight = 400;

f = figure('Position',[0,0,nc * figWidth, nr * figHeight]);
for i = 1:n
    switch i
        case 1 
            color = 'b';
            lsty = '-';
        case 2
            color = green;
            lsty = ':';
        case 3
            color = 'r';
            lsty = '-.';
        case 4
            color = grey;
            lsty = '--';
    end
    if (i == 1)
        subplot(nr,nc,1)
        a = scatter(PredictedDefaultProb, MeanSpread,'o','MarkerEdgeColor', 'k','MarkerFaceColor', white);
        hold on
        plot(PredictedDefaultProb, DataLinearFit,'k-.');
        b = scatter(pDPartition, ModelMeanSpread(:,i), 'd','MarkerEdgeColor',color);
        plot(pDPartition, ModelLinearFit(:,i),'Color',color);
        xlim(xl)
        title('(a) baseline vs. data')
        ylabel('average spread')
        xlabel('default probability')
        hleg = legend([a,b],...
            ['data (int = ',num2str(round(intercept,1)),', slope = ',num2str(round(slope,2)),')'],...
            ['model (int = ', num2str(round(ModelIntercept(i),1)),', slope = ',num2str(round(ModelSlope(i),2)),')']);
        set(hleg, 'location','best')
    end
    
    subplot(nr,nc,2)
    if (i == 1)
        plot(PredictedDefaultProb, DataLinearFit,'k-.');
    end
    hold on
        
    plot(pDPartition, ModelLinearFit(:,i),'Color',color,'Linestyle',lsty);
    hold on
    
end

xlim(xl)
ylim([10,40])
title('(b) compare model variants')
ylabel('average spread')
xlabel('default probability')
hleg = legend('data','baseline','short term','low int. cost','bankruptcy only');
set(hleg, 'location','northwest')

X = [ModelIntercept(1),ModelSlope(1);...
        ModelIntercept(2),ModelSlope(2);...
        ModelIntercept(3),ModelSlope(3);...
        ModelIntercept(4),ModelSlope(4)];%         ModelIntercept(5),ModelSlope(5)];

saveas(f,strcat(FigName,'.eps'),'epsc')






%% Figure A.1: Spreads by income quartile
XDataFileQ1 = strcat(DataPath,'dataforfigures_Sept2024_incomeq1.csv');
XDataFileQ2 = strcat(DataPath,'dataforfigures_Sept2024_incomeq2.csv');
XDataFileQ3 = strcat(DataPath,'dataforfigures_Sept2024_incomeq3.csv');
XDataFileQ4 = strcat(DataPath,'dataforfigures_Sept2024_incomeq4.csv');
FigName = strcat(FigPath,'Fig_A1_Spreads_income');


% load in data
% 0.5 percentile intervals based on default risk
X = readtable(XDataFile);
X = X(1:end-1,:);           % eliminate the last bucket which is a catch-all outlier
X1 = readtable(XDataFileQ1);
X1 = X1(1:end-1,:);           % eliminate the last bucket which is a catch-all outlier
X2 = readtable(XDataFileQ2);
X2 = X2(1:end-1,:);           % eliminate the last bucket which is a catch-all outlier
X3 = readtable(XDataFileQ3);
X3 = X3(1:end-1,:);           % eliminate the last bucket which is a catch-all outlier
X4 = readtable(XDataFileQ4);
X4 = X4(1:end-1,:);           % eliminate the last bucket which is a catch-all outlier


Bins = X.even_fine;
PredictedDefaultProb1 = 100 * X1.mean_default_predict_2;
MeanSpread1 = 100 *  X1.mean_wgt_spread;    % multiply by 100 to convert to percentage points
PredictedDefaultProb2 = 100 * X2.mean_default_predict_2;
MeanSpread2 = 100 *  X2.mean_wgt_spread;    % multiply by 100 to convert to percentage points
PredictedDefaultProb3 = 100 * X3.mean_default_predict_2;
MeanSpread3 = 100 *  X3.mean_wgt_spread;    % multiply by 100 to convert to percentage points
PredictedDefaultProb4 = 100 * X4.mean_default_predict_2;
MeanSpread4 = 100 *  X4.mean_wgt_spread;    % multiply by 100 to convert to percentage points

% rescale to the right level of default
MaxDefault = 20;
i = gridlookup(length(Bins), Bins, MaxDefault);
Bins = Bins(1:i);
PredictedDefaultProb1 = PredictedDefaultProb1(1:i);
MeanSpread1 = MeanSpread1(1:i);
PredictedDefaultProb2 = PredictedDefaultProb2(1:i);
MeanSpread2 = MeanSpread2(1:i);
PredictedDefaultProb3 = PredictedDefaultProb3(1:i);
MeanSpread3 = MeanSpread3(1:i);
PredictedDefaultProb4 = PredictedDefaultProb4(1:i);
MeanSpread4 = MeanSpread4(1:i);


% FIGURE
f = figure('Position',[0,0,600,400]);
xl = [0,MaxDefault];

plot(PredictedDefaultProb1, MeanSpread1,'k-s')
hold on 
plot(PredictedDefaultProb2, MeanSpread2,'b-.d');
plot(PredictedDefaultProb3, MeanSpread3','r:o');
plot(PredictedDefaultProb4, MeanSpread4,'Color',green,'Linestyle','--','Marker','x');
xlim(xl)
title('spread vs predicted default by income quartiles')
ylabel('average spread (annualized pp)')
xlabel('default probability (\%)')
hleg = legend('1st quartile income', '2nd quartile income','3rd quartile income','4th quartile income');
set(hleg, 'location','best','box','off')

saveas(f,strcat(FigName,'.eps'),'epsc')


%% Figure A.2: Spreads by revolver status
XDataFile = strcat(DataPath,'dataforfigures_Sept2024.csv');
XDataFile1 = strcat(DataPath,'dataforfigures_Sept2024_transactors.csv');
XDataFile2 = strcat(DataPath,'dataforfigures_Sept2024_revolvers.csv');
XDataFile3 = strcat(DataPath,'dataforfigures_Sept2024_revolvers_light.csv');
XDataFile4 = strcat(DataPath,'dataforfigures_Sept2024_revolvers_heavy.csv');
FigName = strcat(FigPath,'Fig_A2_Spreads_RevolverStatus');


% load in data
% 0.5 percentile intervals based on default risk
X = readtable(XDataFile);
X = X(1:end-1,:);           % eliminate the last bucket which is a catch-all outlier
X1 = readtable(XDataFile1);
X1 = X1(1:end-1,:);           % eliminate the last bucket which is a catch-all outlier
X2 = readtable(XDataFile2);
X2 = X2(1:end-1,:);           % eliminate the last bucket which is a catch-all outlier
X3 = readtable(XDataFile3);
X3 = X3(1:end-1,:);           % eliminate the last bucket which is a catch-all outlier
X4 = readtable(XDataFile4);
X4 = X4(1:end-1,:);           % eliminate the last bucket which is a catch-all outlier


Bins = X.even_fine;
ActualDefaultProb = 100 * X.mean_default_actual_2;
PredictedDefaultProb = 100 * X.mean_default_predict_2;
MeanSpread = 100 *  X.mean_spread_weighted;    % multiply by 100 to convert to percentage points
MeanSpread_unweighted = 100 *  X.mean_spread_unweighted;    % multiply by 100 to convert to percentage points

% SDSpread = 100 * X.sd_spread;
PredictedDefaultProb1 = 100 * X1.two_pred_mean;
MeanSpread1_wgt = 100 *  X1.spread_wgt;    % multiply by 100 to convert to percentage points
MeanSpread1 = 100 *  X1.spread_mean;    % multiply by 100 to convert to percentage points; this is unweighted by balances

PredictedDefaultProb2 = 100 * X2.two_pred_mean;
MeanSpread2_wgt = 100 *  X2.spread_wgt;    % multiply by 100 to convert to percentage points
MeanSpread2 = 100 *  X2.spread_mean; 

RevolverShare = X.revolver_share;

% rescale to the right level of default
MaxDefault = 20;
i = gridlookup(length(Bins), Bins, MaxDefault);
% PredictedDefaultProb = PredictedDefaultProb(1:i);
PredictedDefaultProb1 = PredictedDefaultProb1(1:i);
MeanSpread1 = MeanSpread1(1:i);
MeanSpread1_wgt = MeanSpread1_wgt(1:i);
PredictedDefaultProb2 = PredictedDefaultProb2(1:i);
MeanSpread2 = MeanSpread2(1:i);
MeanSpread2_wgt = MeanSpread2_wgt(1:i);

RevolverShare = RevolverShare(1:i);
MeanSpread_all = MeanSpread1 .* (1-RevolverShare) + MeanSpread2 .* RevolverShare;
 
% FIGURE
f = figure('Position',[0,0,600,400]);
xl = [0,MaxDefault];

plot(PredictedDefaultProb2, MeanSpread2_wgt,'k-s')
hold on
plot(PredictedDefaultProb2, MeanSpread2,'b-.d');
plot(PredictedDefaultProb1, MeanSpread1,'r:o');
plot(PredictedDefaultProb1, MeanSpread_all,'Color',green,'Linestyle','--','Marker','x');

xlim(xl)
title('spread vs predicted default by revolver status')
ylabel('average spread (annualized pp)')
xlabel('default probability (\%)')
hleg = legend('revolvers weighted','revolvers unweighted','transactors unweighted','all unweighted');
set(hleg, 'location','best','box','off')

saveas(f,strcat(FigName,'.eps'),'epsc')


%% Figure A.3: Figure 1 for broader default measure 
FigName = strcat(FigPath,'Fig_A3_BroaderDefault');
X = readtable(XDataFile);
X = X(1:end-1,:);           % eliminate the last bucket which is a catch-all outlier


Bins = X.even_fine;
ActualDefaultProb = 100 * X.mean_default_actual_3;
PredictedDefaultProb = 100 * X.mean_default_predict_3;
Count = X.count;
ShareBalances_Revolved = X.share_revolved_balance;
ShareBalances = X.share_balance;
AvgBalance = X.avg_revolved_balance;
ShareRevolvers = 100 * X.revolver_share;
ShareRevolversQ1 = 100 * X.q1_revolver;
ShareRevolversQ2 = 100 * X.q2_revolver;
ShareRevolversQ3 = 100 * X.q3_revolver;
ShareRevolversQ4 = 100 * X.q4_revolver;

% rescale to the right level of default
i = gridlookup(length(PredictedDefaultProb), PredictedDefaultProb, MaxDefault);
Bins = Bins(1:i);
ActualDefaultProb = ActualDefaultProb(1:i);
PredictedDefaultProb = PredictedDefaultProb(1:i);
Count = Count(1:i);
ShareAccounts = 100 * Count / sum(Count);
ShareBalances = 100 * ShareBalances(1:i) / sum(ShareBalances(1:i));
ShareBalances_Revolved = 100 * ShareBalances_Revolved(1:i) / sum(ShareBalances_Revolved(1:i));
AvgBalance = AvgBalance(1:i);
ShareRevolvers = ShareRevolvers(1:i);
ShareRevolversQ1 = ShareRevolversQ1(1:i);
ShareRevolversQ2 = ShareRevolversQ2(1:i);
ShareRevolversQ3 = ShareRevolversQ3(1:i);
ShareRevolversQ4 = ShareRevolversQ4(1:i);

% FIGURE
nr = 1;
nc = 3;
f = figure('Position',[0,0,nc * 450, nr * 300]);
xl = [0,MaxDefault];


subplot(nr,nc,1)
plot(PredictedDefaultProb, ActualDefaultProb,'k-o');
hold on 
plot(PredictedDefaultProb,PredictedDefaultProb,':','Color',grey);
xlim(xl)
title('(a) actual vs predicted default')
ylabel('realized default rate (\%)')
xlabel('predicted default probability (\%)')
hleg = legend('realized default rate',...
    '$45^\circ$ line (predicted = realized)');
set(hleg, 'location','best','box','off')

% panel (b): share of accounts
subplot(nr,nc,2)
plot(PredictedDefaultProb, cumsum(ShareAccounts),'k-o')
hold on
plot(PredictedDefaultProb, cumsum(ShareBalances),'b-s')
plot(PredictedDefaultProb, cumsum(ShareBalances_Revolved),'r-.d')
ylabel('cumulative share (\%)')
xlabel('predicted default probability (\%)')
title('(b) distribution of accounts and balances')
hleg = legend('of accounts','of cycle-end balances','of revolved balances');
set(hleg, 'location','best','box','off')


subplot(nr,nc,3)
plot(PredictedDefaultProb, ShareRevolversQ1,'k-s')
hold on
plot(PredictedDefaultProb, ShareRevolversQ2,'b-.d')
plot(PredictedDefaultProb, ShareRevolversQ3, 'r:o')
plot(PredictedDefaultProb, ShareRevolversQ4,'Color',green,'Linestyle','--','Marker','x');
xlim(xl)
title('(c) share of accounts revolving')
xlabel('predicted default probability (\%)')
ylabel('share (\%)')
hleg = legend('at least 1 qtr', 'at least 2 qtrs','at least 3 qtrs','at least 4 qtrs');
set(hleg, 'location','best','box','off')

saveas(f,strcat(FigName,'.eps'),'epsc')



%% Figure A.4: Figure 2 for broader default measure
FigName = strcat(FigPath,'Fig_A4_Spreads_BroaderDefault');
XDataFile = strcat(DataPath,'stats_even_groups_full_fine_revolved.csv');

X = readtable(XDataFile);
X = X(1:end-1,:);           % eliminate the last bucket which is a catch-all outlier

Bins = X.even_fine;
ActualDefaultProb = 100 * X.mean_default_actual_2;
PredictedDefaultProb = 100 * X.mean_default_predict_2;
ActualDefaultProb_Broad = 100 * X.mean_default_actual_3;
PredictedDefaultProb_Broad = 100 * X.mean_default_predict_3;
MeanSpread = 100 *  X.mean_spread_a;    % multiply by 100 to convert to percentage points
SDSpread = 100 * X.sd_spread;


i = gridlookup(length(Bins), Bins, MaxDefault);
Bins = Bins(1:i);
ActualDefaultProb = ActualDefaultProb(1:i);
PredictedDefaultProb = PredictedDefaultProb(1:i);
ActualDefaultProb_Broad = ActualDefaultProb_Broad(1:i);
PredictedDefaultProb_Broad = PredictedDefaultProb_Broad(1:i);
MeanSpread = MeanSpread(1:i);
SDSpread = SDSpread(1:i);

% linear fit
one = ones(i,1);
mdl = fitlm(PredictedDefaultProb_Broad,MeanSpread);
intercept = mdl.Coefficients.Estimate(1);
intercept_SE = mdl.Coefficients.SE(1);
slope = mdl.Coefficients.Estimate(2);
slope_SE = mdl.Coefficients.SE(2);
AdjustedR2 = mdl.Rsquared.Adjusted;


% FIGURE
f = figure('Position',[0,0,1200,400]);
nr = 1;
nc = 2;
xl = [0,MaxDefault];

subplot(nr, nc, 1)
a = plot(PredictedDefaultProb_Broad, MeanSpread,'k-o');
hold on
b = plot(ActualDefaultProb_Broad, MeanSpread,'b:');
c = plot(PredictedDefaultProb_Broad, intercept + slope * PredictedDefaultProb_Broad,'r-.');%,'Color',col3);
xlim(xl)
text(14,15.5,{...
    ['intercept ',num2str(round(intercept,2)),...
    ' (',num2str(round(intercept_SE,2)),')'],...
    ['slope ',num2str(round(slope,2)),...
    ' (',num2str(round(slope_SE,2)),')']},...
    'Color','red','FontSize',10)

title('(a) broader definition of default')
ylabel('average spread (annualized pp)')
xlabel('default probability (\%)')
hleg = legend([a,b,c],'average spread vs. predicted',...
    'average spread vs. realized',...
    ['linear fit (adj. $R^2$ = ',num2str(round(AdjustedR2,3)),')']);
set(hleg, 'location','southeast','box','off')


subplot(nr, nc,2)
basevalue = min(MeanSpread - SDSpread);
area(PredictedDefaultProb, MeanSpread + SDSpread,'FaceColor','k','EdgeColor','k','FaceAlpha',.1,'Linewidth',0.01,'basevalue',basevalue);
hold on
area(PredictedDefaultProb, MeanSpread - SDSpread,'FaceColor',white,'EdgeColor','k','FaceAlpha',1,'Linewidth',0.01,'basevalue',basevalue)
a = plot(PredictedDefaultProb, MeanSpread,'k-o');
xlim(xl)
yl = ylim;
ylim([basevalue,yl(2)])
title('(b) baseline default measure w/ spread error bands')
ylabel('average spread (annualized pp)')
xlabel('default probability (\%)')
hleg = legend(a,'average spread vs. predicted, $\pm$ 1 SD');
set(hleg, 'location','southeast','box','off')

saveas(f,strcat(FigName,'.eps'),'epsc')


%% Figure A.5: Chargeoff (done elsewhere!)


%% Figure C.1: Decisoins without extreme value shocks
FigName = strcat(FigPath,'Fig_C1_DSSY');
NoEVPath = strcat(BasePath, 'NoEV\');

% load in model 
% baseline
mu0 = f2m(strcat(BasePath,'mu0'), [na, ni, ne1]);
pA0 = f2m(strcat(BasePath,'pA0'), [na, na, ni, ne1]);
pD = f2m(strcat(BasePath,'pD'), [naNeg, ni, ne1]);
pD7 = f2m(strcat(BasePath,'pD7'), [naNeg, ni, ne1]);
pD13 = f2m(strcat(BasePath,'pD13'), [naNeg, ni, ne1]);
pDR = f2m(strcat(BasePath,'pDR'), [naNeg, ni, ne1]);

% no EV 
gai0 = f2m(strcat(NoEVPath,'gai0'), [na, ni, ne1]);
gD = f2m(strcat(NoEVPath,'gD'), [naNeg, ni, ne1]);
gD7 = f2m(strcat(NoEVPath,'gD7'), [naNeg, ni, ne1]);
gD13 = f2m(strcat(NoEVPath,'gD13'), [naNeg, ni, ne1]);
gDR = f2m(strcat(NoEVPath,'gDR'), [naNeg, ni, ne1]);


% analysis
cumShareDefault_Base = zeros(naNeg,1);
cumShareDefault_NoEV = zeros(naNeg,1);

cumShare7_Base = zeros(naNeg,1);
cumShare7_NoEV = zeros(naNeg,1);

cumShare13_Base = zeros(naNeg,1);
cumShare13_NoEV = zeros(naNeg,1);

cumShareR_Base = zeros(naNeg,1);
cumShareR_NoEV = zeros(naNeg,1);

DR_Base = sum(mu0(1:naNeg,:,:).* pD,[1,2,3]);
DR_NoEV = sum(mu0(1:naNeg,:,:).* gD,[1,2,3]);

DR7_Base = sum(mu0(1:naNeg,:,:).* pD .* pD7, [1,2,3]);
DR7_NoEV = sum(mu0(1:naNeg,:,:).* gD .* gD7, [1,2,3]);

DR13_Base = sum(mu0(1:naNeg,:,:).* pD .* pD13, [1,2,3]);
DR13_NoEV = sum(mu0(1:naNeg,:,:).* gD .* gD13, [1,2,3]);

DRR_Base = sum(mu0(1:naNeg,:,:).* pD .* pDR, [1,2,3]);
DRR_NoEV = sum(mu0(1:naNeg,:,:).* gD .* gDR, [1,2,3]);

for ai = 1:naNeg
    xBase = 0;
    xNoEV = 0;
    
    x7Base = 0;
    x7NoEV = 0;
    
    x13Base = 0;
    x13NoEV = 0;
    
    xRBase = 0;
    xRNoEV = 0;
    
    for i = 1:ni
    for e1i = 1:ne1
        xBase = xBase + mu0(ai,i,e1i) * pD(ai,i,e1i) / DR_Base;
        xNoEV = xNoEV + mu0(ai,i,e1i) * gD(ai,i,e1i) / DR_NoEV;
        
        x7Base = xBase + mu0(ai,i,e1i) * pD(ai,i,e1i) * pD7(ai,i,e1i) / DR7_Base;
        x7NoEV = xNoEV + mu0(ai,i,e1i) * gD(ai,i,e1i) * gD7(ai,i,e1i) / DR7_NoEV;
        
        x13Base = xBase + mu0(ai,i,e1i) * pD(ai,i,e1i) * pD13(ai,i,e1i) / DR13_Base;
        x13NoEV = xNoEV + mu0(ai,i,e1i) * gD(ai,i,e1i) * gD13(ai,i,e1i) / DR13_NoEV;
        
        xRBase = xBase + mu0(ai,i,e1i) * pD(ai,i,e1i) * pDR(ai,i,e1i) / DRR_Base;
        xRNoEV = xNoEV + mu0(ai,i,e1i) * gD(ai,i,e1i) * gDR(ai,i,e1i) / DRR_NoEV;
    end
    end
    cumShareDefault_Base(ai) = xBase;
    cumShareDefault_NoEV(ai) = xNoEV;
    
    cumShare7_Base(ai) = x7Base;
    cumShare7_NoEV(ai) = x7NoEV;
    
    cumShare13_Base(ai) = x13Base;
    cumShare13_NoEV(ai) = x13NoEV;
    
    cumShareR_Base(ai) = xRBase;
    cumShareR_NoEV(ai) = xRNoEV;
end

nr = 2;
nc = 3;
f = figure('Position',[0, 0, nc * 500, nr * 300]);
subplot(nr,nc,1);
plot(aGrid(1:naNeg), cumShareDefault_Base)
hold on
plot(aGrid(1:naNeg), cumShareDefault_NoEV)
xlabel('debt $a$')
ylabel('share of all default')
title('(a) all default')
hleg = legend('baseline','no EV');
set(hleg,'location','best')

subplot(nr,nc,2);
plot(aGrid(1:naNeg), cumShare7_Base)
hold on
plot(aGrid(1:naNeg), cumShare7_NoEV)
xlabel('debt $a$')
ylabel('share of Ch 7')
title('(b) Ch 7')

subplot(nr,nc,3);
plot(aGrid(1:naNeg), cumShare13_Base)
hold on
plot(aGrid(1:naNeg), cumShare13_NoEV)
xlabel('debt $a$')
ylabel('share of Ch 13')
title('(c) Ch 13')

subplot(nr,nc,4);
plot(aGrid(1:naNeg), cumShareR_Base)
hold on
plot(aGrid(1:naNeg), cumShareR_NoEV)
xlabel('debt $a$')
ylabel('share of delinquency')
title('(d) delinquency')

x = [DR_Base;DR7_Base;DR13_Base;DRR_Base];
y = [DR_NoEV;DR7_NoEV;DR13_NoEV;DRR_NoEV];
X = [x,y]


% savings
aDecision_Base = zeros(na, 1);
aDecision_NoEV = aDecision_Base;
for ai = 1:na
    apVal_Base = 0;
    apVal_NoEV = 0;
    
    muTot = sum(mu0(ai,:,:),[1,2,3]);
    
    for i = 1:ni
    for e1i = 1:ne1
        muVal = mu0(ai,i,e1i) / muTot;
        
        apVal_Base = apVal_Base + sum(aGrid.* pA0(:,ai,i,e1i)) * muVal;
        apVal_NoEV = apVal_NoEV + aGrid(gai0(ai,i,e1i)) * muVal; 
    end
    end
    
    aDecision_Base(ai) = apVal_Base;
    aDecision_NoEV(ai) = apVal_NoEV;
end

subplot(nr,nc,5);
plot(aGrid, aDecision_Base);
hold on 
plot(aGrid, aDecision_NoEV);
xlim([aGrid(1),5])
xlabel('wealth $a$')
ylabel('average decision $a^\prime$')
title('(e) differences in $a^\prime$ choices')

saveas(f,strcat(FigName,'.eps'),'epsc')



%% Figure C.2: Ex ante vs ex post default probability stuff
FigName = strcat(FigPath,'Fig_C2_ExAnte_ExPost');

% load in model 

nAhead = 4;
nPartition = 40;
pDPartition = linspace(0,20,nPartition+1);
TotalBorrowBin = zeros(nPartition,n);
IntegratedCondDefProbBin = zeros(nPartition,n);

for m = 1:n
    q = f2m(strcat(Paths{m},'q'), [naNeg,na,ni,ne1]);
    pA0 = f2m(strcat(Paths{m},'pA0'), [na, na, ni, ne1]);
    mu0 = f2m(strcat(Paths{m},'mu0'), [na, ni, ne1]);

    ExpDefProb = f2m(strcat(Paths{m},'ExpDefProb'), [na, ni, ne1, nAhead]);
    ExpDefProb = 100 * ExpDefProb(:,:,:,4);

    CondDefProb = f2m(strcat(Paths{m},'Cond_DefProb'), [naNeg, ni, ne1, nAhead]);
    CondDefProb = 100 * CondDefProb(:,:,:,4);

    IntegratedCondDefProb = zeros(na, ni, ne1);
    ProbBorrow = IntegratedCondDefProb;
    TotalBorrow = ProbBorrow;
    for e1i = 1:ne1
    for i = 1:ni
    for ai = 1:na
        ProbBorrow(ai,i,e1i) = sum(pA0(1:naNeg,ai,i,e1i));
        TotalBorrow(ai,i,e1i) = mu0(ai,i,e1i) * dot(pA0(1:naNeg,ai,i,e1i), -aGrid(1:naNeg));
        IntegratedCondDefProb(ai,i,e1i) = ...
            dot(pA0(1:naNeg,ai,i,e1i), CondDefProb(:,i,e1i)) / ProbBorrow(ai,i,e1i);
    end
    end
    end

    for j = 1:nPartition
        for e1i = 1:ne1
        for i = 1:ni
        for ai = 1:na
            if (ExpDefProb(ai,i,e1i) > pDPartition(j) && ExpDefProb(ai,i,e1i) <= pDPartition(j+1))
                TotalBorrowBin(j,m) = TotalBorrowBin(j,m) + TotalBorrow(ai,i,e1i);
                IntegratedCondDefProbBin(j,m) = IntegratedCondDefProbBin(j,m) + IntegratedCondDefProb(ai,i,e1i) * TotalBorrow(ai,i,e1i);
            end
        end
        end
        end
    end
end
IntegratedCondDefProbBin = IntegratedCondDefProbBin ./ TotalBorrowBin;

% linear fits
LinearFit = zeros(nPartition, n);
for j = 1:n
    mdl = fitlm(pDPartition(2:end),IntegratedCondDefProbBin(:,j));
    intercept = mdl.Coefficients.Estimate(1);
    slope = mdl.Coefficients.Estimate(2);
    
    LinearFit(:,j) = intercept + slope * pDPartition(2:end);
end

% create figure
nr = 1;
nc = 1;
f = figure('Position',[0,0,nc*600,nr*400]);
for i = 1:n-1
    switch i
        case 1
            color = [0,0.5,0];
            lsty = ':';
        case 2
            color = 'r';
            lsty = '-.';
        case 3
            color = grey;
            lsty = '--';
    end
    subplot(nr,nc,1)
    if (i == 1)
        plot(pDPartition(2:end),LinearFit(:,1),'b');
        hold on
    end
    plot(pDPartition(2:end),LinearFit(:,1+i),'Color',color,'LineStyle',lsty);
end
plot(pDPartition(2:end),pDPartition(2:end),'k--','Linewidth',1);
scatter(pDPartition(2:end),IntegratedCondDefProbBin(:,1),'d','MarkerEdgeColor','b');
xlabel('ex ante default prob. (\%)')
ylabel('avg. ex post default prob. (\%)')
hleg = legend('baseline','short term','low int. cost','bankruptcy only','$45^\circ$');
set(hleg,'location','best','box','off')

saveas(f,strcat(FigName,'.eps'),'epsc')





%% end of program